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Abstract 

We present the constrained Monte Carlo (CMC) algorithm for the QCD evolu- 
tion. The constraint resides in that the total longitudinal energy of the emissions 
in the MC and in the underlying QCD evolution is predefined (constrained). This 
CMC implements exactly the full DGLAP evolution of the parton distributions in 
the hadron with respect to the logarithm of the energy scale. The algorithm of the 
CMC is referred to as the non-Markovian type. The non-Markovian MC algorithm 
is defined as the one in which the multiplicity of emissions is chosen randomly as the 
first variable and not the last one, as in the Markovian MC algorithms. The former 
case resembles that of the fixed-order matrix element calculations. The CMC al- 
gorithm can serve as an alternative to the so-called backward evolution Markovian 
algorithm of Sjostrand, which is used for modelling the initial-state parton shower 
in modern QCD MC event generators. We test practical feasibility and efficiency 
of our CMC implementation in a series of numerical exercises, comparing its results 
with those from other MC and non-MC programs, in a wide range of Q and x, down 
to the 0.1% precision level. In particular, satisfactory numerical agreement is found 
with the results of the Markovian MC program of our own and the other non-MC 
program. The efficiency of the new constrained MC is found to be quite good. 
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1 Introduction 



It is commonly known that the evolution equations of the quark and gluon distributions 
in the hadron, derived in QCD by using the renormalization group or diagrammatic 
techniques [Ij, can be often interpreted probabilistically as a Markovian process; see for 
instance the review of ref. [2 . 

The Markovian type Monte Carlo (MC) algorithm^, implementing the QCD or QED 
evolution equations, is the basic ingredient in all parton-shower-type MCs. An uncon- 
strained forward Markovian MC, with the evolution kernels from perturbative QCD/QED, 
can only be used for the final-state radiation (FSR). It is dramatically inefficient for the 
initial-state radiation (ISR), because the hard process selects certain values of the parton 
energy and the type of the proton. The MC algorithm that allows us to fix (predefine) the 
energy and the type of the parton exiting the parton shower, entering the hard process, 
we shall call the constrained MC algorithm or simply the constrained MC (CMC). 

For the ISR parton shower an elegant example of the constrained MC is the backward 
Markovian evolution MC algorithm of Sjostrand which is a widely adopted effective 
solution in all popular MC event generators, notably in PYTHIA jH] and HERWIG j^j- 
However, the backward Markovian evolution algorithm is a kind of work-around, because 
it does not solve, strictly speaking, the QCD evolution equations. It merely exploits their 
solutions coming from an external, typically non-MC, program. 

The problem that we are posing and solving in the present work is the following: Is it 
possible to invent an efficient constrained MC algorithm, not necessarily of the Markovian 
type, which solves internally the full QCD DGLAP evolution equations on its own, without 
relying on the external non-MC solutions? Since this is a highly non-trivial technical 
problem in the area of the MC techniques, it is necessary to clearly state the motivations 
that have led us to posing it and investing quite some effort in finding at least one 
satisfactory solution. The most important motivation is that we hope to gain more power 
in the modelling of the ISR parton showers - this can be potentially profitable for the 
better integration of the complete next-to-leading-logarithmic (NLL) corrections in the 
parton shower MCs. We also hope for an easier MC modelling of the unintegrated parton 
distributions Dk{Q,PT,x) and MC modelling of the CCFM-type evolution [8J in QCD. 



Let us define more precisely the terminology to be used in this work, since it varies 
in the literature quite a lot. By Markovian MC algorithm we understand an algorithm in 
which the number of emissions (determining the dimension of the phase-space integral) is 
generated as the last variable or one of the last ones^. On the other hand we shall call a 
non-Markovian MC algorithm the one in which the number of emissions (the dimension 
of the integral), is generated randomly as one of the first variables. 

^In the literature, see 0, the Markovian process is usually understood as an infinite walk in a parameter 
space, in which the consecutive steps, indexed by means of the time variable, obey the rule that every 
single step forward is independent of the past history of the walk. 

^The backward Markovian algorithm is also very well described in ref. p]- 

•^Our definition of the Markovian MC implies that the corresponding Markovian process is terminated 
by means of some stopping rule, for example by limiting the process time or other parameter. 
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The name of "constrained MC" can also be associated with a MC in which the dis- 
tribution to be generated is restricted to a less-dimensional hyperspace by inserting the 
6{F{xi, ...,Xn)) function, where F(xi,...,x„) = defines the constraint. The CMC al- 
gorithm should efficiently generate points within this subspace. It is widely known that 
it is usually much more complicated to generate the distribution on such a hyperspace 
that in the entire space, even if the original unconstrained distribution is simple, for ex- 
ample it is the product of many simple distributions. The energy-momentum conserving 
§W(^P — '^pi) is a well known example of such a constraint. 

In our case the constraint is given hy S{x — Y[zi), where x is predefined (the outermost 
integration variable, generated in a CMC as the first one) and guments of the 

DGLAP evolution kernels. We shall describe an example of the CMC algorithm which 
fulfills such a constraint, measure its efficiency and make a numerical test of its correctness, 
by comparing its results with those of the traditional unconstrained Markovian MC, and 
other non-MC programs. 

The outline of the paper is the following: in section 2 we introduce the basic notation 
and rederive an iterative solution of the DGLAP evolution equations. Since the full so- 
lution is algebraically involved, we discuss in section 3 our CMC solution for the simpler 
case of pure bremsstrahlung out of quark or gluon. Section 4 describes the hierarchical 
reorganization of the iterative solution, which is necessary for the full scale solution with 
arbitrary number of gluon-quark transitions. This complete solution, along with numer- 
ical tests, is presented in section 5. Summary and outlook concludes the paper. For the 
sake of completeness, a small appendix lists the standard LL QCD kernels. 

This paper describes the main part of a wider effort on the MC modelling of the 
QCD evolution equations. In ref. |9j the reader may find an alternative CMC algorithm 
for the QCD evolution, which looks slightly inferior to the one presented here and is 
implemented numerically only for pure bremsstrahlung. References [10^ and ^T] are 
devoted to a precision MC modelling of the LL and NLL DGLAP evolution equations 
using the Markovian (unconstrained) class of algorithms. In this context it is worth to 
remind the reader that the consistent integration of the NLL perturbative corrections 
at the fully exclusive level in the MC event generator of the parton shower type is the 
challenging problem both theoretically and practically. For the recent efforts in this 
direction see refs. [T2t IT8] , for example. The algebraic proof of certain important identity 
used in this work is published separately in ref. [Hj. 

2 Iterative solution of the QCD evolution equations 

Let us rederive the iterative solution of the QCD evolution equations. The starting point 
is the DGLAP |lj set of evolution equations: 

1 

^^D,{t,x) = J2 I ^P>.,{z)^D,[t,^)=J2^,,{t,.)^D,{t,.), (1) 
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where 



dxidx26{x - xiX2)f{xi)g{x2), 



(2) 



and "ykjit, z) = ^^^Pkj{z). Indices i and k = G,q,q denote gluon or quark - the evolution 
time is t = ln{Q). The differential evolution equation can be turned into the following 
integral equation: 



3*fc{t,t0) 



t 

Dk{t,x)=Dk{to,x) + J c?tie*'=(*i'*°)^?|(ti,-)®l^,(ti,-)(a^) 



(3) 



to 



where e is the IR regulator in the kernels, 



'?%{t,z) = 7k,{t,z)Q{l-z-e), 
and the Sudakov form factor is given by 

Mt,to)= [ dt'7i,{t',e). 

J to 



(4) 
(5) 



(6) 



The complete set of LL kernels (Pfcj is collected in an appendix. 

The multiple iteration of the above integral equation leads us to an iterative solution of 
the evolution equations, given by the following series of integrals, ready for the evaluation 
with the help of the standard MC methods: 



t 1 

oo r- n p p 

Dk{t, x) = e-*'=(*'*°)Z}fc(to, n y - ^-i) J ^-^^ 

n=l fco...fc„_i '- i=l Q 



(7) 



i=l 



Dko{to,Xo)6lx ~ XoY[z, 



i=l 



where fc„ = k. In ref. ^U] the above equation was solved with the three-digit precision 
using the Markovian MC method^. In this work the above solution is the starting point 
for constructing the CMC algorithm. 

There is still one standard technical point: the one- loop dependence of the strong 
coupling as(t) = 27r/(/5o(t — InAo)) may destroy the efficiency of the MC algorithm, 
unless it is conveniently compensated by means of the mapping tj — Tj = ln(tj — InAg). 

^The energy sum rules J2i I ^-^ z'-^ik{z) = are instrumental in this MC method and in fact the 
equation for xDk{x) rather than D}~{x) is solved there. 
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Once it is done, the iterative solution transforms into the following form: 



°o pi r " pT pL 

n=l ko...kn-i '-*=! "^^0 



X e 



-*fe(-r,r„) 



■*fci_l(T!,Ti-l) 



1- i=l 



-Dfco(ro,xo)5( X - xo 



i=l 



(8) 



where k = kn- The kernel 7 and form factor $ are from now on redefined as follows: 

as{t) dt 2 

-'kiki-ii^i, ^i) = ~ ^kiki^ii^i) = ^ ^kiki^xK^i) ^ 



$fc(r,ro) 



(9) 



TO 



In the present LL case, J* becomes completely independent of Tj. See also refs. 13 IH] for 
more discussion on the optimal choice of the evolution time. 

In the following we shall often employ the short-hand notation: 



6x>y = 9{x — y), where 9x>y = 1 for x > y, otherwise O^yy = 0, 
6x=y instead of 6{x — y), 



(10) 



so as to keep formulas more compact. 

Throughout this work we concentrate on the LL case. However, our CMC algorithm 
solves most of the problems on the way to the CMC solution for the NLL DGLAP. The 
exact Monte Carlo solution of the NLL DGLAP evolution equations in the framework of 
the unconstrained Markovian approach is presented in a separate work (ref. ^Ij). It may 
serve as a numerical benchmark for the future work on the NLL CMC. 



3 CMC for pure bremsstrahlung only 

Before we unfold all the details of our construction of the CMC algorithm for the full 
DGLAP equations, let us first describe it for the simpler case of the pure QCD brems- 
strahlung. This will help the reader to follow all algebraic technicalities of the full CMC 
solution. It will also provide a building block for the full CMC solution. 

The starting point is the iterative solution of the QCD evolution equations of eq. (jHl) 
truncated to the multiple gluon bremsstrahlung emitted from the parton k = G,q, q, with 
the starting distribution Dk{To,Xo) = 6{xo — 1): 

x2)fcfc(r,ro;x)=e-*'=(^'^«)|5,=i + f^^ f[ f dr, f dz, zA{t„ Zi)Sx=Yi^^A ■ (11) 
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For the purpose of the MC generation, we introduce simphfied kernels: 

1 , l^ _ ^ ejl -e-z) 

1- z 

^ (12) 



zP^aiz) ^ zP^M = zBGcOil - ^ - ^) ( + i ) = ^GG- 



e/ X . Ae. ^ _ d <d{l - e - z) 



A. A/ 

Note that in eq. ()12p we do not modify the virtual parts of the kernels. We also do not 
assume any sum rules relating J"^^ to / zTlf,. This is so because we view eq. (fTT|) as a 
part of a bigger framework (eq. (jH)) for example) with both quarks and gluons, and we 
will assume later on a complete sum rule "Plf. = / z'J'^j,^. 

The above simplification of eq. ()12j) will be countered by the MC weight wp defined 
below: 

^ ^ ln{l-a;) r 

t „=i i=i , -/ ^ ^ J 



ln(e) 



(13) 



where yi = ln{l — Zi) and 



^kir, To) = (r - To) ( fefc In - - ttfc 



(14) 



_ 2 _ 2 

= -^Bkki flfc = -T^kk- 

Po Po 



The factor x'^'^ is introduced in order to obtain as uniform shape of the MC weight 
distribution for all x G (0, 1) as possible. It will also help to correct for the fact that 
Pqq{zi) does not have a l/zj component - in this case Hi = ^ can be pulled out of the 
integrand. For the moment we assume ujg = and u;^ = 1, but obviously they are dummy 
parameters, which may influence the MC efficiency but not the final MC distributions. 
The energy constraint 



x = X\ Ziiyi) = n ~ exp(|/i)) = F(yi, y2,...,yn), = F{y) or 

1=1 i=l 

1 

- = ^fiVj) = -lni^(y); /(%■) = - In (1 -exp (yj)) = - In Zj 



In 

X . 



provides also an integration limit z{yi) > x, which translates into yi = ln(l — Zi) < 
ln(l — x). Taking advantage of the symmetry of the integrand we can introduce ordering 
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in the y-variables: 

oo n 

where |/o = |/min = ln(e) and |/max = ln(l - x). 

The function is very steeply (exponentially) rising, hence the constraint x = 

Y[i=i ^iiVi) is effectively "saturated" by a single Zj while the other ones are — 1, i ^ j- 
In terms of y variables, i/j ~ ?/max; while the other ones i/i, i ^ j, move freely within 
the (i/min, 2/max) interval. In our case, because we ordered y-variables, it is the biggest 
Un — Umax that effectively takes responsibility for satisfying the constraint. 

In the following we translate the above statements into rigorous mathematics. The 
procedure of eliminating the energy constraint goes in three steps: 

Step one. We introduce the new integration variable Y. Its introduction is immediately 
countered by the 5-function: 



f oo 
^ n=l 



(17) 



Step two. We perform the following simple linear transformation 

y^ = v\ - Y- (18) 

Note that Y was "adjusted" from the very beginning such that = l/n + ^ = l/max, see 
also a graphical illustration in fig. [TJ 



y y i y n y 

min ymax 



Figure 1: Graphical representation of the linear transformation yi y[ 
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The Jacobian of the above hnear transformation is equal to 1. We obtain: 

^ 71=1 



(19) 



where ?/o = ymm + Y, and Oy^^y^.^^ = Oy[>Y+y^in = ^y[>y'o defines the IR lower boundary of 
the phase space. 

Step three. The energy constraint 5{x — F) is eliminated by means of the F-integration: 



n=l 1=1, 



dn wp > , 



(20) 



TO 



where Yq = YqIxjJ') is the solution of the transcendental equation x = F{y' — Y). One 
subtle point is that thanks to y^ = yjam + Yoix, y'), and Yq >0 we may keep formally the 
same integration limits y'^ G (?/min, Z/max) as before. 

Summarizing the above three steps: we effectively traded complicated 6{x — F{y)) into 
much simpler 6{y'^ — |/max)- We may also say that we projected, by means of the parallel 
shift^ (along the diagonal), points from the hyperspace yn = ymax into the hyperspace 
x = F{y). 

Let us now have a closer look into formula for c}yF(y' — Y), which is necessary for the 
MC weight and for the numerical solution of the equation x = F{y' — Y): 



\dYlnF{y'-Y)\Y=Yo 

" I- Zk 



E 

k=l 



J2dym-Y) 

k=l 



Y=Yo 



k=l 



exp(i/fc - Y) 



Zk 



Y=Yo 



(21) 



|Z = 2fc 



k=l 



The additional transformation y'- = y^^^ + r^Ay with A^, = |/max — 2/min and < < 1 
as well as exporting part of the integrand into MC weight brings the integral closer to the 



'^The projection of a similar type (with dilatation instead of parallel shift) has been employed in the 
CMC already in refs. JJ5 , while the idea of saturating the constraint with a simpler one using single 
variable can be traced back to the MC program of ref. (and its unpublished versions) as well as to 
ref. [HI- 
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form suited for the MC generation: 



1 T 

+ x^"-' Yl ^fe^r' n / ^n>n.Arn - 1) / dnw* I, 

n=l i=l r( _ ^ 



TO 

# 



(22) 



The average or maximum weight is improved if the weight wf(zi, Zn) is rescaled by 
the function g{x) (which can be pulled out of the integral). The weight distribution is 
conveniently stabilized if we take for g{x) the biggest term in \dylnz{y)\z=Zk, with the 
largest y^, which we approximate ets Zk — x: 

^* = '^P I a 1^ TP^^^^^'vw ^y[-Yo>ymir.^ 9ix) = \dy\n z{y)\z=a: = (23) 



\dYlnF{y'-Y)\y=Y, 



X 



Finally we introduce the normalized variables Sj = (tj — Tq) / (t — Tq) , rescale w* and 
symmetrize over rf 

a;Dfefc(r,ro;a;) = e-*'=(-'^°)|5,=i+ 

1 1 

n-i j-i Q Q 

A = &fc(r - ro)Aj, = 6fc(r - ro)[ln(l - a;) - Ine] = Di(l - x) - ^{e), 
CR(a;) = 6a;(t — To) In a;. 

Introducing explicitly a Poisson distribution P{n\X) = e~^\^ /n\^ and remembering that 
$fc(T, To) = — (6fe Ine + ak){T — tq), we obtain the following master formula on which the 
MC algorithm is built: 

xD,,(r,ro;a;) = V e(---°)«'= e^(^)5„=o5.=i + Wi-.>ee^(^-^)^^^(r - Tq) 

X p(n - -X)- m) n /^n " / ^''^ 





Tj = To + Si(T - To), 2;^ = 1 - e^* = 1 - exp(ymin + rAy - Yq). 

In our CMC algorithm for pure bremsstrahlung one generates first then n, next 
Zi (constructed from Ui) and finally Tj (constructed from Sj). The algorithm is clearly 
non-Markovian, as is seen from the fact that the number of emissions is generated as the 
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second variable, not as the last one. The distribution of the variable x is done according 
to the following primary distribution: 



I xg{x) J 



(26) 



which is obtained by means of neglecting and performing all summations and inte- 
grations. The MC weight is, of course, restored later on. 

It is quite convenient, in the construction of the MC program, that in the above 
distribution we are able to extend artificially the integration above x = 1 — e, by means of 
mapping x into the new variable U = exp(CR(l — x)), so that we can generate x as if there 
was no 5x=i- This resembles quite strongly the analogous trick for the multiple photon 
emissions in the YFS-type MCs for QED ^H]. Let us work out the details, restoring the 
initial parton distribution at t = and the hard process function H{x): 



1 

dx H{x) Dkk{T, x) 



i 1 

dx j dZ J dxo 5{x — XqZ )H{x) T>'^^{t,tq] Z) Dfc(ro,Xo) 



j dxH{x) I ^Z--V— ")'^45z=i e''^^^ + e,.z>se''^'-'^\dzm-Z)\}D,{ro, 



1 1 

X 

Z 

ei X 

1 exp(K(l-x)) 

JdxH{x) J rft/Z(f/r-V---°)'^'={5^=oe^(^)+ WpWe))}/^fe(ro, 



Z{U) 

€1 

1 exp(3?(l-a;)) 

X 



ZiU) 

£1 

U{Z) = e^(^-^) = (1 - zf^^^-^°\ 
Z{U) = 1 - exp {{hu{T - ro))-^ ln[/), 

(27) 

and remembering that 1 > Z{U) > 1 — e is mapped exactly into one point ai U = 0, 
reproducing the component ~ 6z=i, i.e. /;^P(^(^))t/f/ = exp(a^(5)). Once X is chosen, n is 
generated according to a Poisson distribution (shifted by 1) and then all variables and 
Zi are generated. 

The above collection of detailed formulas determines uniquely the whole CMC algo- 
rithm for the pure bremsstrahlung case^. For the sake of completeness let us summarize 

^The only element that is not described in fine detail is the method of solving the transcendental 
equation for the constraint lo(j/i)- We use a variant of the standard method of the tangentials. It is in 
principle quite straightforward - the only complication is that it must work for all values of y[. For certain 
values, because of dyF ^ 0, attention has to be paid that the number of the iterations is sufficient. 
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point by point the complete CMC algorithm: 

• The outmost integration variable generated as the first one is total x, the argument 
of the hard process cross section H{x). 

• The second generated variable is Z, the total loss of energy due to multiple gluon 
bremsstrahlung, with the help of the mapping: U{Z) = e^(^~^) = (1 — zY''^'^~'^"\ 

• The generation of variables x and U (and of k = k^, if necessary) is done with the 
help of the general-purpose MC tool FOAM [TTJll^ . 

• Knowing Z{U), if Z > 1 — e, the emission multiplicity n is generated according to 
a Poisson distribution Pn~i (non-Markovian!), otherwise Z = 1 and n = 0. 

• Variables Si, i = 1, 2, n are generated uniformly and mapped onto rj(sj) and ti(rj). 
They are ordered. 

• Not ordered variables G (0, 1) are generated, such that one of them is set equal 
to^ 1; they are mapped into y'iiri). 

• The solution Y = Yq of the transcendental equation lnF(y' — F) — Inx = is 
found numerically (NB: the derivative (9ylnF for the MC weight is obtained as a 
byproduct). 

• With Yq at hand, all variables Zi{yi{y'j)), i = 1,2, ...,n are calculated. 

• In the case yi < ymm, see fig. El the MC weight is zero and the algorithm stops. 

• Finally the MC weight lu* is calculated. Optionally the weighted MC event is 
transformed into an unweighted one with the help of the standard rejection method. 




Figure 2: Graphical representation of the rescaling procedure yi — > y'^ 

In this way we completed the description of the CMC algorithm for the pure brems- 
strahlung case, which will be the essential building block for the general-case CMC algo- 
rithm described in the following. 
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Figure 3: Comparison of CMC and EvolMC in the case of pure bremsstrahlung out of gluon. 
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3.1 Numerical test of the CMC for pure bremsstrahlung 

The CMC for pure gluon bremsstrahlung was tested separately for emission from gluon 
and quark lines by comparing its results with those from the Markovian MC EvolMC of 
refs. PO] and ^T]. As in ref. results are for the LL DGLAP evolution of the parton 
distribution in the proton from energy scale Q = 1 GeV to Q = 1 TeV. We use the same 
starting parton distributions in proton at Q = 1 GeV as in ref. [10] and the same range 
of X G (0.001, 1). In the upper plot of fig. El we show the gluon distribution evolved up 
to scale Q = I TeV (lower curve) due to pure gluon bremsstrahlung, obtained using our 
new CMC and the Markovian EvolMC. Results are indistinguishable; we therefore plot 
their ratio in the lower plot of this figure. The results of the two programs agree perfectly, 
within the statistical errors, which are of order 0.5% at low x. MC results are for statistics 
of a few hundred millions of MC events. A higher-statistics comparison will be presented 
in the following section. We also include, in the upper plot of fig. El the result of the 
evolution (from the Markovian MC) due to all transitions, not only bremsstrahlung. As 
we see, the complete result differs from the pure bremsstrahlung one by a factor of almost 
2, hence the inclusion of the transitions of gluon into quark and back is very important. 

In fig.Elwe present numerical results of the analogous comparisons for the quark singlet 
Q = q + q. Again, pure bremsstrahlung results of the new CMC agree very well with 
these of the Markovian EvolMC, to within of the statistical error, which is about 0.5% in 
most of the x range. Contrary to the previous case the curve for full evolution coincides 
with that for pure bremsstrahlung. This is easy to explain as a result of the suppression 
of the gluon distribution at high x, resulting in the smallness of the Q G contribution. 

On the technical side, let us remark that there are two methods of obtaining pure 
bremsstrahlung contributions from the Markovian MC. One may simulate full DGLAP 
evolution, including Q ^ G transitions and select events (evolution histories) in which 
only pure bremsstrahlung occurs. The other method is to suppress kernels for Q ^ G 
transitions completely. In the present version of the Markovian EvolMC, both methods 
are available, and both give identical results. In the present work we mainly use the first 
method of selecting evolution histories out of complete evolution. 

4 The CMC with the flavour transitions - fuU DGLAP 
4.1 General discussion 

Before getting into details, let us describe the essential ingredients of our CMC algorithm 
for full LL DGLAP, with an arbitrary number of flavour-changing transitions G ^ Q, 
where Q = q,q. The first ingredient is the observation, made in ref. JU] , that the average 
number (n) of G ^ Q transitions is much lower than the average number (A^) of G ^ G 
or Q — > Q ones, that is of the gluon bremsstrahlung emissions. In fact (n) ~ 1 for 
the evolution from Qo = 1 GeV to Q = 1 TeV, while (A^) ~ 20 (for e = 10~^). This 
suggests quite strongly that we should consider the evolution process (emission chain) 

''We generate n of them and rescale such that the biggest is equal 1. 
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as a two-level process: sublevel of the pure bremsstralilung and superlevel of the flavour 
transitions. Let us call it "hierarchical" organization of the emission chain. The hope 
is that the small number of superlevel transitions can be modelled, for example by a 
general-purpose MC tool such as FOAM, thanks to a relatively small dimensionality of the 
problem. The assumption is that pure bremsstrahlung segments can be treated separately 
and efficiently. 

As discussed in ref . , for the treatment of the energy constraint (5-function, there are 
at least two options. We may deal with it independently and separately for each of the 
two levels; this is called type I solution in ref. jH]. In the other CMC solution, nicknamed 
CMC type II in ref. [9 , the energy constraint is implemented globally, for both levels at 
once, using the assumption (corrected later on by the MC weight) that Dk{tQ,x) ~ x^*""^. 
The latter solution is described and implemented in ref. [H], up to bremsstrahlung level, 
with the explicit algebraic layout for the full DGLAP. 

In the following we shall walk along the path of CMC class I solution, the superlevel 
implementation employing the general-purpose MC tool FOAM and the sublevel being 
implemented exactly as in the previous section. Both levels feature the energy constraint, 
that is the total loss due to multiple emissions/transitions is predefined and in the MC 
it is generated as one of the first variables. The same is true for the total number n of 
fiavour transitions and the number of gluon emissions n^, in every i-th pure bremsstrahlung 
segment of the emission chain. 

4.2 Hierarchical organization of the emission/evolution chain 

As we have argued in the above discussion, the two-level "hierarchical" reorganization of 
the emission/evolution chain is mandatory for any reasonable CMC scenario. In algebraic 
language, the transition to hierarchical organization means that the sums over the fiavour 
indices in the iterative solution of the evolution equations of eq. (jH)) are reorganized 
in such a way that all adjacent gluon emission vertices are lumped together into the 
distributions/integrals Dkk{T,To', x) of the previous section. The remaining integrals and 
sums will belong to the superlevel. The formal derivation of the resulting hierarchical 
iterative solution of the evolution equation is presented separately in ref. [Hj. In the 
following we shall present only the final result, discuss its structure and apply it to the 
CMC algorithm. 

The iterative solution of the DGLAP equation reorganized in the hierarchical form 
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T > T n k=k 



D (T,X)= 



X =x =,x^Z riz. Z. 

" n+l n+r J J 

X. k. ^ k ^ V X k 



'A / \ ^.^n^ k^^ j^ N VlVl 



1111-1 ' ' r n 



Figure 5: The scheme of kinematics and flavour indices in the hierarchical Markovian. 



reads as follows: 



1 1 



n=l fe„_i...,fci,fco 

fcj7^fej_l,j = l,....n 



TO 



Y\_ j dzi j dZi j dxo 

^—1 n n n 







X 



n 



n 



>- 1=1 



(n n+l \ 
X - Xii^Zi^ZiY 
i=l i=l ^ 

T 1 

n=l 1=1"'"/ 



(28) 

where we denote k = kn and the virtual form factors are defined as follows: 

1 

(29) 

-Rjfc = dz z'J'fi^{z), B!f. = ^ = Rk — Rkk- 
J^'^ 

The (ifc-function obeys the normalization condition Jq dZ Zdk{T, Z\tq) = 1. It is easy 
to check that d^ is related to the pure bremsstrahlug distribution Ti^k worked out in 
section IHl eq. (fTT| . in the following way: 

2)fcfc(r, To; Z) = e-(^--°)«'^4(r, Z\ro), 
Rk = 6fcln(l/£:) - flfc- 



(30) 
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As is also schematically shown in fig. El the integrand of eq. consists of a product 
of the superlevel transition probabilities, each of them consisting of the fiavour- changing 
kernel and the pure bremsstrahlung function. The whole emission chain is terminated 
with yet another bremsstrahlung function. 

Let us add the following interesting side remark: We could interpret the above formula 
as the Markovian process, each step being an entire (pure bremsstrahlung) Markovian 
process of its own. Of course, it would be possible to implement such a hierarchical two- 
level Markovian scenario as the MC simulation of the unconstrained QCD evolution*; 
however, in such a case, contrary to the CMC case, there is no convincing reason to invest 
an extra effort to implement in practice such a solution. 



4.3 The CMC solution type I 



Exploiting the two- level hierarchical solution of eq. ()28|) . we shall work out in the following 
the CMC solution of type I, concentrating on the fiavour-changing superlevel, because the 
pure bremsstrahlung level is the same as described in section 01 except that it is repeated 
here not once but many times in a single evolution process. As a first step, let us write 
eq. ()28j) once again, in a more compact form, eliminating J dxo at the expense of the 
(5-f unction: 



Dk(T,x) 



Xo 



dZ Dkk{r, To; Z) —Dk{To, Xq) 9x<z+ 



X 



+ E E U 



n=l fc„_i...,fei,fco i=l 



1 



X 



dZn+1 Dkk{T, Tn; Zn+l) 



Tj_i; Zi) 



(31) 



X 



Xo = X 



f Zn+1 Y\ ^i^i ] 
^ i=l ' 



Next, we work out the explicit kinematic limits in terms of x- variables defined as follows: 



Xi 



'\^^j^j]xO — xi^n+l ^j^j] ) ^ — 0) 1) 2, .., n, Xn+l = X, Zn+i = 1. 

' j=l ' ^ j=j+l ' 

(32) 



^We tried to check in the hterature whether this possibihty was noticed by the authors of the classic 
Markovian MCs, but we could find no explicit reference to such a scenario. 
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In this way, we obtain 



^ / dZi T'fcfc(r, To; Zi) Xo£'fc(^o, Xq) + 

J X 

T 1 

oo p n „ ~\ P 

n=l fc„_i...,fci,fco '-J=l^„ -' „ 



+1 '^kki'T, Tn] Zn+l] 



kjjtkj_i,j = l,...,n 
1 



" TO 
1 



(33) 



Til Ti-l'i Zi) 



The functions 2)fc._^fc._j(rj, rj_i; Zj) are the multidimensional integrals of their own, de- 
scribed in section |3J which in the Monte Carlo are implemented as a separate module 
providing pure bremsstrahlung subevents with the weight In the first stage of the 

MC, neglecting w = Ilj^C^i' i-^- replacing D^k by 2)^;. of eq. (j26|) . we have to generate 
3n + 1 continuous variables explicitly present in eq. according to the integrand 



(l~x)''k(^—^0) 



Dk(T,x) = X ^ 



dUi z([/i)"^-2e'^'=(^-^°) xoDk{To,xo) + 

(l-x)''k(^-^ri) 



+ 5Z n / d-r^ ^-.>-.-i / dUn+l Z{U, 

n=l fc„_i...,fci,fco i=l_„ n 



kjT^kj_i,j — l,....n 
1 



i=i 

X Xoi5fco(^0,^o), 

where t/j = exp(CR(Zj)) and 



(l-Xi/Zi) 



/ 



-1 



(34) 



(35) 



The first term in = 0) in the above sum is identical to the one discussed in section IHl 
The second term for = 1 is the new and non-trivial one, representing one q G or 
G q transition accompanied by the two segments of the pure bremsstrahlung. It reads 
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as follows: 



T 



(l_a;)6fc(T-Ti) 




-1 



fen ^ 



/ 







(36) 



1 



'^kJn—^o) 



I 



dUi Z{Ui) 



"^fen-SpafeoCn-To) 



xo-Dfco(^o,a;o) 







where xi = X/Z2, xq = xj {Z2Z\Z\) and k\ = k. 

In addition to continuous variables, we have to solve the problem of generating effi- 
ciently all (izscreie variables fcj, z = 0, 1, . . . , n—\. For gluon and 2nf quarks and antiquarks 
in the contribution with n flavour transitions, we have in eqs. ()31|) - (jMj) up to (2?t,j + 1)" 
terms. This might be a serious problem in the case of implementing all these component 
integrals one by one in the general-purpose MC tool such as Foam, taking for example 
^ < ''^max = 4. Luckily, thanks to symmetries valid in the LL approximation, we are able 
to get this problem under control, at least for massless identical quarks. This method- 
ology is not the same as the traditional splitting of parton distributions into singlet and 
non-singlet parts, although it exploits the same properties of the kernels. 

The essence of the solution of the above problem can be demonstrated in a trans- 
parent way for just one type of quark and antiquark. The extension to the case of n/ 
identical massless quarks is not difficult and will be done later on. Let us analyse the 
flavour sum '^k^_-^kn-2 fci fco' ^^^^ condition ki 7^ fcj-i, i = 1,2, ... ,n, in eq. (jMj) . We 
split the parton distribution of eq. (jH^ into components with a well deflned number of 
flavour transitions D(x) = '^^='0 and we shall discuss the components D*-"-* one by 

one. In the above and in the following we omit all bremsstrahlung contributions in D^'^\ 
integrations over Zi, etc., in eq. (j34j) - we keep track of only the flavour- changing kernels 
and their indices. 

The easiest case is n=0, that is the case of the pure bremsstrahlung. Here, the inte- 
grand contains just one term proportional to Dk, modulo the bremsstrahlung part, where 
k is one of the three possible flavours k = G,q, q. In the distribution given to Foam this 
contribution is treated separately and is generated automatically by Foam with the correct 
probability. 

The flrst non-trivial case is that of n = 1 and k = ki = G. Here, the sum under 
consideration is D^"* = Xlfeo -^Gfco-Dfco, where fco 7^ G, hence ko = q,q. Since Tcq = ^Gg 
we may replace both of them by Tcq and pull them out of the sum. We obtain Dq^ = 
Gq J2 ko=q g^f'o — "^GqDq, where we have introduced the inclusive (singlet) initial quark 
distribution Dq = Y.k,=q,qDk,- 

In the similar case k = ki = q for n = 1, the sum under consideration is Dq^'' = 
'^kfj'^qkoDko- In this case the only possible contribution is for k^ = G, and the only 
remaining term is D^^^ = "PqcDG- In the case of tagged antiquark, k^ = q, we obtain 
exactly the same contribution: Df^ = TgcDc. 
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The general case of the n > 1 transitions is quite similar to the n = 1 case. In the 
case of A; = A;„ = G, in the sum 

we may repeat the previous reasoning we followed for n = 1. In the first step wc reduce 
the summation over the last index kn-i = q,(l with the help of "J^cq = '^Gq and = 
obtaining: 

= 27Gq ^ '^qkr,-2 ■ ■ ■ '^k^koDko- 
kn-2---kl,ko 

In the next step we get rid of the sum over kn-2 

Dq^ = ^'PQq'S'qG ^ ^Gfen-a • • • '^kikoDko- 

fc„_3...fcl,fco 

We continue with the elimination of the sums one by one, obtaining for odd n: 

^(n) ^ ^(n-i)/2 y^^y^^y^^ . . . j)^^ Bq, 

while for even n we obtain: 

Df^T^I-" T>Gq'^qGT'Gq...'^qG ^G- 

As a result of the above reasoning, the whole sum is reduced to just one term for each n. 

(Every n = 0, 1, . . . , Wmax is generated by Foam separately.) It is important to stress that, 
in the MC, where we are interested in a fully exclusive history of the intermediate states 
in the emission chain, we may easily undo the summation over equal contributions from 
quarks and antiquarks, leading to factors 2^""^^/^ or 2"/^, and choose randomly with the 
probability 1/2 between kj — q and kj — q for every intermediate non-gluon state (every 
other link in the emission tree). 

The case oi k = kn = q and n > 1 can be analysed in an analogous way: 

'^,k^-.---'^k,k,Dk,. 

kn-lkn-2---kl,ko 

Due to kernel properties and kn ^ kn-i, the condition Dq^^ reduces in the first step to: 

^f^J'^G Yl '^Gkr.-2---'^k^koDko. 

k„-2...ki,ko 

The final result is either 

^(n) ^ ^(n-m y^^y^^y^^ . . . jJq, for n odd 

or 

^(n) ^ 2"/2 TqC-^Gq^qG ■■■'^PqG Dq, for U CVCU. 
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Summarizing, for identical quarks, in the LL approximation, we may effectively get 
rid of summations over all of the flavour indices! The case of n/ identical quarks is quite 
analogous - the only difference is that we get (2n/)^"~-^^/^ or (2nj)"/^ weight factors in 
front of each single final term. When restoring the type of the intermediate quark we 
choose randomly with equal probability one of the 2nf quarks and antiquarks. 

The above collective treatment of the intermediate quarks/antiquarks in the process 
of implementing the distributions of eq. (jH^ as integrands of Foam is relatively easy for 
the finite number of flavour transitions (n = 0, 1, 2, 3, 4). As we remember, the only case 
where we need to explicitly generate an individual quark or antiquark type ko according 
to ~ Tfco(xo)i^fco, is the case of n = 0, i.e. pure bremsstrahlung. In all other cases, 
n > 0, we may treat quarks and antiquarks at the intermediate stage of the MC gener- 
ation collectively and randomly choose their individual type (index) later on, with equal 
probabilities. 

In the above explicit algebra, the difference between the traditional split into singlet 
and non-singlet components D'^ = Dq + Dq and D'^^ = Dq — Dg and our alternative 
split into pure bremsstrahlung component D^^ or D'^^ and the rest Dg^^'^^ = JJ^^^^^ is 
manifest. Both techniques exploit, of course, the same symmetry properties of the kernels. 




-3 -2.5 -2 -1.5 -1 -0.5 

log (X) 




log (x) 



Figure 6: CMC versus Markovian MC for gluons; number of quark-gluon transitions n = 
0, 1, 2, 3, 4 and the total. The ratio in the lower plot is for n = 0, 1 and the total (blue). 
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Figure 7: CMC versus Markovian MC for quarks; number of quark-gluon transitions n 
0, 1, 2, 3, 4 and the total. The ratio in the lower plot is for n = 0,1 and the total (blue). 



With all the above algebra at hand, we can formulate our CMC algorithm in the 
general LL DGLAP case: 

• Generate superlevel variables n, ki, Ti Zi and Zi using the Foam general-purpose MC 
tool according to eq. 

• In the above we limit the number of flavour transitions (G — > Q and Q ^ G) to 
n = 0, 1, 2, 3, 4, aiming at a precision of ~ 0.2%. 

• For each pure gluon bremsstrahlung segment defined by Zi and (rj,rj_i), i = 
1,2, ...,n + 1, gluon emission variables (zj*'', rj*^), j = 1,2, ...,n^*\ are generated 
using the previously described dedicated CMC. 

• Weighted events are generated. They are optionally turned into weight-1 events 
using the conventional rejection method. 

The above algorithm is already implemented in the C++ programming language and 
tested using Markovian MC EvolMC of refs. [T0| and pHj; see next section. In the program 
construction the central variable is the number of quark-gluon transitions n. We do not 
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know the analytical expression for the distribution of this variable. It is one of the variables 
managed by Foam. In the initialization phase, Foam finds out the distribution of n numer- 
ically and then generates n efficiently in the range < n < nmax- In the CMC program 
construction the rimax = 0, 1 cases were programmed first and the results were compared 
with EvolMC. The main programming exercise consist of the careful programming of the 
integrand distribution for Foam. It was finally done for arbitrary nmax- The other part 
of the programming is related to joining pieces of several pure bremsstrahlung segments 
into a single long emission chain with an arbitrary number of quark-gluon transitions. 
Object-oriented programming tools made this task easier. 



4.4 The CMC for DGLAP: numerical results 

The basic tests of the new CMC algorithm are presented in figs.lHlandlZl As in section ITTI 
we examine results of the DGLAP evolution from the energy scale Q = 1 GeV to Q = 1 
TeV, using as the starting point quark and gluon distributions in proton dX Q = 1 GeV, 
exactly the same as in ref. PUI- In fig- El we show distributions of gluon at Q = 1 TeV, 
while in fig. [3 are plotted the results for quark, Q = g + g, at the same high scale Q = 1 
TeV. In these two figures we compare gluon and quark distributions obtained from the 
new CMC program and from the Markovian EvolMC^ of ref. jTUj. The main numerical 
results from both programs, marked as "total" in the upper plot of both figures, are 
indistinguishable. We therefore plot their ratio in the lower plot of both figures. They 
agree perfectly well within the statistical error, in the entire range of x. For x < 0.1 the 
statistical error is below 0.1%. 

The plots in figs. El and [71 contain, however, more tests than that for the total nor- 
malization. As already mentioned, in the process of constructing the CMC program we 
have tested also each "slice" of the gluon and quark distribution for a given number of 
quark-gluon transitions n = 0, 1, 2, 3, 4, on the way from 1 GeV to 1 TeV. For example, 
in fig. El we show separately the contributions to the gluon distribution from the following 
evolution histories: 
G 

G and any number of gluon emissions out of Q and G, 
Q G, etc. 
G^Q^G, etc. 

Q ^ G ^ Q ^ G, etc. ("Total" is the sum of n = 0,1,2,3,4.) 
They are shown in the upper plot of the figure one by one for the two programs compared, 
CMC and EvolMC. As before, the results are indistinguishable - this is why we plot their 
ratios in the lower plot of the figure. The discrepancy is within the statistical error, as 
in the total contribution. In this plot the comparison of the two programs for the n = 
slice is a repetition of the pure bremsstrahlung test from section 13.11 but for much higher 
statistics. 



n 


= 


G 


n 


= 1 


Q 


n 


= 2 


G 


n 


= 3 


Q 


n 


= 4 


G 



^This Markovian MC program has been tested in refs. ^01 and against two non-MC programs 
QCDnumie [2] and APCheb33 j22j. Small systematic discrepancy between EvolMC and QCDnuml6 for gluon 
distribution reported in ref. JOl was later eliminated and explained in ref. |11|. 
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Analogous slices for the quark distribution: 
n = 0: Q^Q 

n = 1: G Q and any number of gluon emissions out of Q and G, 
n = 2: G ^ Q ^ G ^ Q, etc. 
n = 2: G ^G ^Q, etc. 

n = 4: Q^G^Q^G^Q, etc. ("Total" is the sum of n = 0, 1, 2, 3, 4) are shown in 
fig- HI Again the results from the new CMC and Markovian EvolMC agree perfectly well. 




Figure 8: Weight distribution of the new CMC algorithm for gluon (upper plots) and quark 
(lower plots) as a function of a;. 

In fig. ISl we show the distributions of the MC weight distribution separately for the 
evolution yielding gluon and quark at Q = 1 TeV. This is an important technical test of 
the MC integration/simulation performed with the help of the general-purpose MC tool 
Foam. As we see, the resulting MC weight is well limited, w < 1, and the average weight 
is about 0.08. The efficiency is therefore quite good, substantially better than that for 
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method II. b of ref. jH]- Note that for nmax = 4 the integrand of Foam is 15- dimensional, 
but the modelling of this distribution is performed using only 2000 cells. This clearly 
demonstrates the power and the usefulness of this tool. 

5 Conclusions and outlook 

In this work we have shown that there exists an efficient constrained Monte Carlo algo- 
rithm for the initial-state radiation emission chain in QCD. Its most likely application 
will be in the construction of a new class of parton shower Monte Carlo event generators 
for the QCD initial-state radiation. 

In this context it is worthwhile to mention that a similar CMC algorithm of class I 
has been worked out [SB] for the HERWIG-type QCD evolution, see ref. i.e. for the 
z-dependent — z)Q) and t-dependent IR cut-off e(t)^°. 

The present CMC solution is restricted to the LL kernels. It will be interesting to ex- 
tend it to the NLL case. The preparatory step in this direction is already done. In ref. ^1] 
the MS NLL kernels are implemented within the Markovian Monte Carlo program. They 
can be ported to the non-Markovian CMC in the future, if necessary. 

Let us stress that the CMC algorithm of this work is not the only one known. For 
instance an alternative non-Markovian CMC algorithm class II exists; see refs. and 
[0]. It is defined there for the full DGLAP, although it is implemented/tested for the pure 
bremsstrahlung only. It has worse MC efficiency than the algorithm presented here and 
leads to higher dimensionality of the integrand managed by FOAM. 

Another possible application of the presented CMC algorithm is the MC modelling 
of the unintegrated parton distributions, including these based on the CCFM-type evo- 
lution^^, for the purpose of MC simulating W and Z production at LHC. In fact, the 
unintegrated parton distributions Dk{kT,x) are already calculated |26J from the one-loop 
type CCFM model of ref. j2Z], in the framework of unconstrained Markovian MC. 
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^°This has been done for the case of pure gluonstrahlung combined with any number of quark-gluon 
transitions. For more technical details see also http: / /home. cern.ch/jadach/public/desy _mar05.pdf 
^^For example in the approach similar to that in the CASCADE Monte Carlo of ref. |25) . 
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Appendix: QCD LL kernels 

We present here a table of the elements in the LL kernels {Tf = njTfi), Q = q + q 



IK 


AO) 


r(0) 


^(0) 
^IK 




Dik{z) 


JdzDf^iz) 


GG 


11 2rTl 


2Ca 


2Ca 


2Ca{-2 + z-z^) 
2Tf{z^ + (1 - zf) 





11^ 


QG 









2Tf 




QQ 




2Cf 





Cf{-1-z) 





-|Cf 


GQ 






2Cf 


Cf{-2 + z) 








Pik{z) = 5(1 - z)6ikAkk + 



— ^^ikBkk H — Cik 
- z 



Dik(z). 



Temporary simplifications for the purpose of the MC generation are: 



5, 



z>e 



GG- 



zPf^{z) ^ zPg{z) 



B, 



'l-Z>£ 



zTtk^t. z) - z) = ^ zP^,{z) 



2B 



kk 



01- 



■z>e 



TC 



/?0(t-tA) l-Z 



'?ik{t) = ^lBkkln--Akk 



(37) 



(38) 
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